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ABSTRACT 


A technique for the use of semi-infinite planes to 
approximate the log magnitude versus log frequency charac- 
teristics of two-dimensional filters is presented (Extended 
Bode approach). This technique is applied to quarter plane 
filters and works well for separable transfer functions. 
Other non-separable canonic (basic) transfer functions are 
also studied in terms of their planar approximation. The 
general techniaue is shown to be useful for the insight it 
provides as well as being a simple approach to design. The 
double bilinear z-transform is studied for use with the 
two-dimensional analog transfer functions to convert them 


to the digital domain. 
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I. INTRODUCTION 


A. DIGITAL SIGNAL PROCESSING 

caca ls Tgnhal processing has its roots in 16th century 
mathematics, especially in the fields of astronomy and the 
compilation of mathematical tables. Today it has become 
a powerful tool in a multitude of diverse fields of science 
and technology. The applications of digital signal pro- 
cessing varies from low-frequency spectrum seismology 
through spectral analysis of speech and sonar into the 
video spectrum of radar systems [1]. 

In their book, Theory and Applications of Digital 
Signal Processing [2], Rabiner and Gold have combined digi- 
tal signal processing theory, with a variety of applications 
ranging from sonar, radar, communication, music, seismic 
and medical signal processing, with digital component 
technology. This technology is the main driving force for 
progress in this field as well as the general area of com- 
puter design. It was also observed by them, that, although 
the formulation of engineering problems is often as vague 
as those of the "softer" sciences (such as anthropology, 
psychology, etc.), the execution of these problems appears 
to depend on greater and greater accuracy and reproducibility. 
The capability of digital systems to achieve a guaranteed 
accuracy and essentially perfect reproducibility is very 


appealing to engineers. 
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A digital filter is defined in [3] to be a computa- 
tional process in which a sequence of numbers acting as 
input, is transformed into a second sequence of numbers, 
or output digital signal; where the term "digital" implies 
that both time (the independent variable) and amplitude are 
quantized. 

The field of one-dimensional digital filtering, which 
is outlined in figure 1.1, encompasses recursive, non- 
recursive and Fast-Fourier Transform processing. The terms 
E urvessnd nonrecursive, instead of IIR (Infinite Impulse 
eec ponse) and FIR (Finite Impulse Response) are preferred 
in this thesis. It is shown that recursive processing is 
much more efficient than nonrecursive processing.  Stockham's 
method [4] to perform fast convolution, which later became 
known as the FFT method, improved the efficiency of non- 
recursive techniques, so that comparisons in one-dimension 
are no longer strongly biased toward recursive methods [2]. 

One-dimensional recursive digital filter design depends 
strongly on the effective and well developed continuous 
filter design theory. Moreover, the stability analysis of 
higher order one-dimensional recursive realizations can be 
solved by investigation of the root distribution of its 
factored form representation. It is important to realize 
that continuous domain design techniques and factorization 
property (the fundamental theorem of algebra) exists only 


in one-dimension. 
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Areas of application of one-dimensional digital signal 
processing are listed in figure l.l. There are important 
military and specifically naval applications in addition 
to those listed, such as missile and torpedo guidance, 
launcher control system, fire control, combat information 
collection, dissemination and display. 

Digital filtering in several dimensions, which is over- 
viewed in figure 1.2, has gained considerable importance, 
especially for the two-dimensional case. Figure 1.3 presents 
an example of two-dimensional low-pass and high-pass filtering. 
Until 1966, two-dimensional digital filtering was implemented 
by nonrecursive or convolutional techniques. In this method 
the output is the weighted sum of unit sample responses of 
all past input values. The serious disadvantage of the 
convolutional method is the requirements of a very large 
number of arithmetic operations. 

The development of the Fast Fourier Transform (FFT) in 
1966, reduced the number of arithmetic operations consider- 
ably and is used extensively today. Filtering via FFT is 
Eleconplished by computing the transform of the input func- 
ren, multiplying by the frequency response of the filter, 
and inverse transforming the result. The recursive 
algorithm, which has in general an infinite impulse response 
[3], constitutes another technique for the realization of 
two-dimensional digital filters. 

Hall [5] compares the amount of computation required EQU 
the three filtering techniques. He shows that, in general, 


the FFT and recursive algorithms are preferable to the 
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An example of two-dimensional filtering 
(An original photograph) [2] 
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Fig. 1.3.b. An example of two-dimensional filtering 
(Low pass filtered version of Fig. 1.3a) 


[2] 
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Fig. 1.3c. An example of two-dimensional filtering 
(High pass filtered version of Fig. 1.3a) 
[2] 


a 





nonrecursive one in number of computations and storage 
requirements. Hall demonstrates that the recursive filter 
algorithm constitutes the best method for large data, i.e., 
it is the fastest and cheapest. 

When processing two-dimensional data by recursive filters 
a fundamental problem exists due to inherent feedback, namely, 
the problem of numerical stability. Since, in general, two 
variable polynomials cannot be factored into a product of 
first and second order real coefficient polynomials in each 
of the variables, it is difficult to solve the stability 
problem for two-dimensional recursive filtering. Conse- 
quently, the majority of papers published in two-dimensional 
digital filtering deal with the design of nonrecursive 
filters which are inherently stable [6-10]. 

There are several papers discussing two-dimensional 
recursive digital filters. For example [11] and [12] 
formulate two-dimensional recursive filters by the z- 
transform and linear difference equations, and although 
they investigate problem areas related to stability and 
realization, the majority of problems remain to be solved. 

The most important applications of two-dimensional 
digital filtering of digital data are in picture processing 
and geophysical data analysis. Picture processing can be 
categorized into: 

a) Digital image restoration and enhancement, 


b) Computer pictorial pattern recognition. 
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The methods of image restoration and enhancement are 
applied to invert degradations, such as aberrations, 
atmospheric effect, scanning, motion, and to manipulate 
images to improve viewing phenomena experienced by the human 
eye. Image restoration and enhancement has been used with 
great success in biomedicine, i.e., in extraction of quan- 
EBBENEIve information from x-ray films, chromosome counting, 
measuring the extent of arteriosclerosis from arteriograms, 
Io rensic sciences application, i.e., fingerprint image 
ennancement for automatic classification, and in astronomy, 
i.e., removal of turbulence from astronomical photography 
a3]. 

Pictoral pattern recognition by digital computer has 
gained great importance in satellite surveillance and mili- 
tary reconnaissance. The application of two-dimensional 
digital filters vermits the separation of different horizontal 
scans on magnetic and gravity maps in geophysics and can 
be used equally well for structural and topographic maps 
or for any other type of factor which is available in the 
Brmecesot a Planar grid [14]. Although digital filtering 
in several dimensions is gaining importance in medicine, 
i.e., heart volume measurements, in fire control problens, 
and to analyze complex electronic circuitry, there exists 
no general theory concerning structures, analysıs and 
design. 

These considerations can be summarized as follows. 


Since one-dimensional recursive digital analysis and design 
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techniques depend strongly on the fundamental theorem of 

wo Dra and ón extensive utilization of continuous design 
Ebeory, it E a Special, i.e., a non-generalizable 
field in the area of N-dimensional recursive digital filtering. 
Multi-dimensional digital filtering is well developed in two- 
dimensions but is based, due to the absence of recursive theory 
analysis and design tools, predominantly on nonrecursive 
filtering schemes. Because of their inherent advantages, much 
of the current research is directed towards recursive filtering 


algorithms and techniques. 


B. AREA OF INVESTIGATION AND OBJECTIVES 

One area of investigation of this thesis is to see how the 
semi-infinite straight line approximation technique of one- 
dimensional filters (based upon the log modulus/log frequency 
or decibel vs. log-frequency plots) can be extended to two 
(or higher) dimensional recursive filters. The advantages of 
log-modulus approach is that a transfer function in factored 
form becomes a linear sum of terms when the logarithm of its 
magnitude is taken, and each of these terms can be approxi- 
mated by semi-infinite straight line segments in the one- 
dimensional case, and semi-infinite planes in the multi- 
dimensional case, when expressed in terms of log-frequency. 

The following are specific objectives: 

Ios cu late typical (canonic) factors for 

two-dimensional filters and to investigate their properties 


in the log-modulus/log frequency domain. 
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Zep lO develop an approximation technique for the 
design of two-dimensional recursive filters to meet given 
Specifications in the frequency domain. 

The frequency response of two dimensional filters is 
inherently a volume in the space of magnitude, log 03 and 
log "Pu This volume is the given or desired specification 
of the filter, which in the separable case can be approxi- 
mated by the combination of semi-infinite planes. 

In this thesis we investigate how we can approximate 
this volume by semi-infinite planes as an extension of 
the approximation used in one dimension, i.e., by semi- 
infinite lines. After determining the planar equations 
which perform the given specification of the filter, we 
investigate ways to get a transfer function; first, in the 
continuous frequency domain and second for two-dimensional 
recursive equations in the discrete variable domain. 

The results Of this study will classify those charac- 
teristics that can be achieved by cascading simple canonic 
factors, and indicate limitations that are imposed by 
trying to achieve a particular frequency characteristic 
(such as the fan filter) in as simple form as possible. 

The approach taken here is to investigate the problem 
Ene eontinuous domain first, and then study the trans- 
formation from continuous frequency (differential equation) 


domain to discrete time (sampled data) domain. 
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SEN PREEVEEN OF RESULTS 

We have proposed a geometric method to design two- 
dimensional recursive digital filters. This is a new 
method which is relatively simple and easy with respect to 
other design methods. We have investigated properties of 
some typical (canonic) factors and considered their stability 
properties. In Chapter II we go through the basic properties 
of two-dimensional filtering. In Chapter III we provide 
background in terms of planar geometry for the discussion 
in subsequent chapters. In Chapter IV we investigate the 
properties ome ical factors and their stability in continu- 
ous variable domain. The new design technique, including 
stability considerations, starting with the continuous 
frequency domain and extending to two-dimensional recursive 
filters in the discrete domain, is presented in Chapter V. 
The results are summarized and some comments for further 
study are made in Chapter VI. It should be noted that the 
results presented apply to quarter plane filters, as defined 


in Chapter V. 
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IE FUNDAMENTALS OF TWO-DIMENSIONAL RECURSIVE FILTERING 


De INTRODUCTION 

There are many signals that are inherently two-dimensional 
for which two-dimensional signal processing techniques are 
required. The most common examples are, of course, images 
in which the two variables are the spatial coordinates. 
Image processing [15], [16] plays an important role in many 
areas of scientific and technical research. Some examples 
are satellite imagery radiography, radar, and biomedical 
images such as acoustical holograms, medical x-rays, and 
electron micrographs [2]. Not all two-dimensional data 
come from an image. In prospecting for gas and oil, a 
common procedure is to monitor the seismic signals produced 
by detonating an explosive charge [17]. An array of geo- 
phones situated along a line radiating in a vertical or 
horizontal direction from the point of explosion is used 
to record these signals. In this case, the two variables 
are distance (from the point of explosion measured along 
the line) and time (the time at which the seismic signal 
reaches a given geophone or, equivalently, a given distance 
from the explosion). 

Although two-dimensional signals may be processed by 
one-dimensional systems, it is generally preferable to con- 
sider using two-dimensional systems.‘ Many of the basic 


ideas of one-dimensional signal processing may be readily 
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extended to the two-dimensional case. There are some 

very important concepts of one-dimensional systems, however, 
that are not directly extentable to two dimensions. It 

is the goal of this chapter to discuss the basic ideas and 
techniques of two-dimensional signal processing and to 
illustrate them into the context of two-dimensional filter 
design. A more detailed discussion of this theory plus 

an overview of recent work in two-dimensional filtering is 
contained in the excellent survey paper by Mersereau and 
Dudgeon [18], and H. Chang and J.K. Aggarwal [19]. Other 
useful references are the first few chapters of [l1] and 


Siapter 7 of [2]. 


D DEFINITIONS 
Ix Dwo-Dimensronal Signals 

One-dimensional signals are functions of a single 
variable and two-dimensional signals are functions of two 
integer variables. Consider the two-dimensional sequence 
x(m,n) where m and n are integer variables, shown in figure 
2.1. AS in the one-dimensional case, the notation x(m,n) 
is often short hand for a sampled version of a continuous 


two-dimensional signal x(s,t): i.e., 


us) = x (mT,,nT,) = x(s,t) s=nT, ‚t=nT, (2. 1) 


One interpretation of x is a function which assigns a 


(generally complex) number to each integer ordered pair 
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(m,n). The function x is not defined when either or both 
ers arguments (m,n) is not an integer. Such a signal 
will be interchangeably referred to as either a two-dimensional 
array, a two-dimensional sequence or simply as a two-dimensional 
signal. 

Some useful two-dimensional sequences are defined 
below and shown in figures 2.2 to 2.4. These include: 


ie Digital Impulse or Unit Sample 


2 0 elsewhere 


il m,n > 0 
(m,n) = (23) 
0 elsewhere 


m n 
a] a» mue O 
x(m,n) = (2.4) 
0 otherwise 
4. Sinusoid (Complex) 
j (o4mton) 
x(m,n) = e —o « m,n € +œ (2.5) 


As seen above, the two-dimensional step is related to the 


two-dimensional impulse by 
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m 
u_, (m,n) = dl ) uj Uns ,ni) (2.5a) 


2.  Two-Dimensional Systems 
A two-dimensional discrete system is defined as 
the unique transformation T that maps an input sequence 
{x(m,n)} and the initial condition sequence {s(m,n) } into 


an output sequence {y(m,n)} 

y(m,n) = T[{x(m,n)},is(m,n) }] (2 10D 
In signal processing applications, we assume zero initial 
conditions. With this assumption, the above relation may 
be shortened to 


une = Tlixim,n);] o 


The system characterized by T is 


(a) linear if 
T[fax, (m,n) + bx, (m,n) }] = aT({x, (m,n) }] + bT[{x, (m,n) 3] 
(2.8) 
for arbitrary constants a and b. 
e shift invariant if 
y(mm,,n-n,) = T[(x (m-m, ,n7ni) 7] o) 


for all my and m 
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A hardware and software system which is equivalent to the 
transformation T of (2.6) is called a digital filter. The 
class of complex-valued sequences which are square summable 
will play the role of input and output signals of the 
digital filter.  Linearity and shift-invariance are indepen- 
dent properties of a system in that neither property implies 


the other. Thus the system 

T(x(m,n)] = him,n) x(m,n) 
is linear but not shift invariant and the system 

2 

T{x(m,n)] = x (m,n) 

Me shift invariant but not linear 
For two-dimensional discrete linear systems, the 

impulse response is defined as the response of the system 
at (m,n) to a two-dimensional digital impulse input at 
(m4 , ni) DB zero initial conditions, 1.8, 

h(m,n;m, ,ny) = T[(u, (m-m, ,n7ni) j] (2.10) 
where u, (m,n) denotes the digital impulse. In particular, 


the impulse response of a two-dimensional discrete LSI 


Memear Shift Invariant) system is 
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h (m,n;m, ,ni) = h(m-m, ‚n-n,) 


Or 


mane) T{tu, (m,n) }] (2240) 


The response of a two-dimensional LSI*system with zero 
initial conditions is completely characterized by its 


Inoulse response. In fact using the identity 


oo 


x(m,n) = ) Xx (m, ,n,) u5 (mcm, ,n7nj) 
zu PA 
= ) x (m-m, ,n-n4)ug (m4 ,n4) 
m,;,n,--9? 
we get 
Ann) = Te) x (m, ,n4) ug (mm, ,n-n4) )] 
m,,n47-7 
= ) x (m, ,n4) T[ (up (m-m, ,n-n3) 1] 
m4,n4-77 
E ) x (m, ,nj)h(m-m,,n-nj) 
m,n = 
“757 = Linear shift Invariance 
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oo 


y(m,n) = > x (m-m, ‚n-n,)h(m, ‚n,) 


n E 
= A e (m, n) PRR E) 
where "*" denotes two-dimensional convolution. In other words, 


for LSI systems the basic convolution theorem is valid. 
x 
NESST, SCpatabilitcy, Stability 


In discussing various features of two-dimensional 
discrete systems, 1t is convenient to have spectral nota- 
tion [19] to represent certain regions of the spatial domain. 


l. Let S(a,8) denote a sector with the angular 


interval (a,8). In polar form S(a,8) may be described by 
IO, APO y (213) 
where (a,8) are measured in radians. In particular, we 


give the following notation for the sectors frequently 


encountered: 
See SA0, n2]; S-+ = S[t/2,7] 
S-- = S[-1,-1/2]; S+- = S[-1/2,0] (2.14) 
S*+  S[0,17); EE Zen) 

pere [and "]" imply including the boundary and ")" implies 


mote including the boundary. 


ab 





E ceothat $T, S-t, S--, and S+- represent first, second, 
third, and fourth quadrants, respectively, and S*+ and Sx*- 
represent half-planes which are symmetric to each other with 
respect to the origin. (See fig. 2.5.) 

2. For a two-dimensional sequence {x(m,n)}, the 


set 


E, = { (m,n) 





Zr mr 05 (2.15) 


is called the support of {x(m,n)}. 

3. Let {h(m,n) } be the impulse response of a two- 
dimensional discrete LSI system, and let E, be the support 
EXCUh(m,n)]. 

Then, a two-dimensional LSI system is said to be 
causal ifE, is the set S++, and semi-causal if E, is the 
EEUCUSxT. They are shown in Fig. 2.5. 

A two-dimensional discrete LSI system is said 


to be separable if its impulse response can be factored into 


a product of one-dimensional responses (Fig. 2.6); i.e., 


Am,n) ^= h4 (m) h, (n) (2.16) 


Tf the equation (2.16) is not satisfied, the filter is said 
to be nonseparable. The advantages of separable filters is 
that two dimensional convolution (2.12) may be carried out 
as a sequence of one-dimensional convolutions. This can be 


seen by rewriting (2.16) as follows. 
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ig. 2.6. Separable Digital Systems 
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I 
a 


ymn) by (m)h) (n,)x (m-m, ,n-n3) 2,173) 


= ) o hj(m)f 2  hz(ny)x(mm,n-n,)] (2.17b) 


E ) hi (m, ) a (m-m, ,n) (2 STO) 


where a (m-m, ,n) is a sequence of one-dimensional convolutions 
(obtained by evaluating the terms inside the bracket of 
mar’) for each fixed value of m). Equation (2.17c) shows 
that y(m,n) may be obtained by a second sequence one-dimensional 
senvolution. 

If both the input sequence x(m,n) and the filter 
impulse response h(m,n) are separable, then it is readily 
seen that the output sequence y(m,n) is also separable. 


In this case we obtain the result [using (2.12) and (2.16) ] 


y(m,n) = [ } hy, (m,)b(m-m,)][ ) A aa Ie) 
m., ==% n +72 
= a(m)8 (n) (2.18b) 
where 
ma) = bin ein) 


D> 





A two-dimensional LSI system is said to be (BIBO) 
stable if every bounded input sequence produces a bounded 
output sequence. This condition is satisfied if and only 
if the two-dimensional impulse response of the system is 


absolutely summable, i.e., 


) Jn(m,n)| « e (2.19) 


m,n=-» 


As in the one-dimensional case, (2.19) can be shown to be a = 
necessary and sufficient condition for stability [20]. One 
Meoblem with (2.19) is that it can be quite difficult to 
evaluate for an arbitrary h(m,n). 
4. Two-dimensional Difference Equations 

As in the one-dimensional case, two-dimensional LSI 
systems can often be described by a two-dimensional linear 
constant-coefficient difference equation relating the output 
Mn) ot the filter to its input x(m,n). The most general 


form for such a difference equation is shown by 


i alk]),K,)x(m-k),n-k,) - ) b(t k3) y (m-2, ,n-2,) 


(Kk) /k5) eR, (2, 25) ER, 
(2.20) 


where RA and Ry are finite sets of spatial grid points, 


called the input and output mask, respectively, and {a(k,,k,) } 


and {b(2,,25) 3 are the set of constant coefficients that 
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characterize the particular filter. Generally, a difference 
equation has a family of solutions as is the case with a 
differential equation. The difference equation is said to 
be recursive if there exists a sequence of computations which 
yields the output sequence serially from the input sequence 
and that portion of the output sequence which is already 
eenpureqd, including the initial condition. An implicit 
EXdering of grid points in the spatial domain is assumed for 
the serial computation of the output sequence. 

There are two cases of (2.20): 

li Rov Ro in the S++ (first quadrant), b(0,0) # 0 
EM n9 9 for outside of S++ (zero initial conditions), 
then (2.20) represents a causal system. 

DLE Ros Ry in Sx+t, b(0,0) z 0 and y(m,n) = 0 for 
Eu-gce' cf S«4- (zero initial conditions), then (2.20) 
represents a semi-causal system. 
In both cases we may assume b(0,0) = 1 without loss of 


generality and write (2.20) as 


y(m,n) = N a(k, /K,)x(m=k, n-k) 
(Kk) ,K5) eR, 


- ) b(2 
(2, 12)) ER, (0,0) 


(A5) Y (m-£ ;n-25) 


il l 
er 2) 


It is clear that (2.21) is a recursive equation, and the 


implicit ordering, for example may be given as 
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MOON (1,0), E us (2,0) 2 2S XA A 2 32v] 


for a causal system, and 


E779, 1)1,(0,2),...; (sor qom p ys aou... xu cce 


ee core) a (l,2),e. ss ey 


for a semi-causal system. 
IE {b(2,,2,)} = tuo (21,25) ] the equation (2.21) 
becomes 


Amen = i a(K,,kj) x(m-k, ,n-k 2877) 


2! 


so that the output sequence is a weighted moving average of 
the input sequence. The impulse response of (2.22) is 


unique and is obviously {a(k k4)]. Such filters are called 


E 
finite extent impulse filters (FIR) because the impulse 


response has only finitely many nonzero terms. FIR filters 


obviously are always BIBO stable because a(k,,k,) is always 
absolutely summable. (It has only finitely many nonzero 
terms.) 


If, on the other hand, the sequence (b(1,,24)) 
has several nonzero terms, then the impulse response will 
generally have infinitely many terms. We will say that such 


filters are infinite-extent impulse response filters (IIR). 
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Such filters may be BIBO unstable, and so it is necessary 
pssdevise a stability test for such filters. 

FIR filters have the advantages that they are 
easier to design, that controlling the phase response is 
easier (in particular, it is easy to produce filters which 
have a linear phase response), and that since the input- 
output relationship is a finite-extent convolution, Fast 
Fourier Transform (FFT) techniquescan be used to speed the 
computation of the output sequence. On the other hand, 
recursive (IIR) filters generally require less high speed 
storage and usually require less computation time than FIR 
filters. Comparison of FIR and IIR filters is contained 
En [5]. 

5.  Two-dimensional z-tranform 

A useful mathematical tool for representation of a 
Sequence x(m,n) is its two-dimensional z-transform. Given 
a sequence x(m,n), the z-transform of this sequence is 


defined to be: 


A A ID 
Z[ix(m,on)])] 2 X(z,,22,) = D ) x(m,n)z,z, (2.23) 
m=-=w n=-» 
where Z4 and z, are independent complex variables. The most 


Nito roperty Of the z-transform is the fact that the z- 
transform of the convolution of two sequences is the product 


o their z-transforms, i.e.: 
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mirum n)hl*ih(mm] = X (21,25) H(z1 , Z5) (2.24) 


Berefore, the input-output relationship of an LSI filter 


may be expressed in the z-transform domain as: 
Y(z],2,) = H(Z],2,)X(z,,2,) 0225) 


where H(z1:25), the z-transform of the impulse response 
meecalied the transfer function. 

A two-dimensional transfer function evaluated on 
the unit circle is called the frequency response of the 


system, i.e., 


jo jo ju jo jo jo 
Ten een. Ze Le ^) (2.26) 


when a two-dimensional LSI system is describable by the 
difference equation of (2.20), and we take a z-transform 


SEIS, it follows that 


; 1 2 
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Or 


H(z 


1127) Y(2z,),2,)/X(z,),2,) 


) alknık,)z) 2, 
(Ki; K4)ER. 
) BL NN 


= A(z]),2,)/B(z,),2,) 12. 2m) 


where A and B are the z-transforms of ta(k4 ,k5)] and 
(b(24,125)7 respectively. In particular when (2.27) repre- 
Sents a causal or semi-causal system, i.e., (2.27) is the 


ELEranstorm of (2.21), we can give explicit descriptions for 


Ra and Ry in each case. For the causal case we define 
Meg. 2.5a) 
R, = {(k,-k,)|0<k, <M,, 0<k,<M)} 
(23278) 
Ry = {((2,,2,)|0<52, < Ny, 0<2,<N)) 


Ber che semi-causal case we define (Fig. 2.5b) 


} 


aJ 
| 


{(k,,k,) }O<k) <M,,k, -0;-M, Sk £M2,0<k,<M, 


(229) 


U 
il 


( (2412) |O« 244 «Ng, 45 7 0;-N. € 44 <N O< ta SNo} 
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and in both cases b(0,0) = 1 is assumed. The inverse 


z-transform is defined as 


Amn) ë = I f e: nge 


nee 
(2n3)% 
dh 2 


Zo dz,dz. 2230) 


where Cy and cj are Suitable closed contours in the Z4 and 
Zo planes. 

The two-dimensional Fourier transform of a sequence 
is defined as 


Jw Jw 
X (e we 5 wees (cme) | 


=e CZ ae Ze 
W? 2 
(Z4 rZ) ET 


7 -j (o, m*on) 
- )  x(m,n)e (2432) 


my n=-ow 


For any given sequence the set of (24,25) for which che 
z-transform converges is called the region of convergence. 
Uniform convergence of the z-transform requires that the 


sequence be absolutely summable, 1.e. 


A A AS (2.32a) 


m,n=-» 
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jui ju 


Ewslarly, the Fourier transform X(e ¡e ) converges 
En rormlyv to a continuous function of Wy and Wo qf 
co 
) Sn E co, 232) 


m,n=-w 


A power series of the form of (2.23) is known as a two- 
dimensional Laurent series in the theory of complex functions 
in two variables. 
6. Two-dimensional DFT 
Let {x(m,n)} be a finite support sequence of the 
Size MxN. Then the two-dimensional discrete Fourier 


transform (DFT) of {x(m,n)} is defined as 


m -j2r (Ek, + 2k.) 
X(k,,k,) = } ) compe 2733) 
m=0 n=0 
sor 0 < Ky il and 0 < Ko « N-1. Similarly, the inverse 


discrete Fourier transform (IDFT) is given by 


45 





X(m,n) = 57 I ME E b à (994) 
2 


Note that the DFT corresponds to sampling the Fourier trans- 


Exc cot MxN points, i.e., 


IEn ale 
X (K4 , K^) = x(e ¡e ) ki K, 
043277 ,03=213 (2575) 
joi Ju, 
More generally, if X(e ¡e ) is the Fourier transform 


of {x(m,n)} whose size is greater than MxN, then the 
MEET of X (ki ,k5)] will result in aliasing. In other words, 
the IDFT of (X (ks ,k,)] will yield one period of {x(m,n)} 


where 


x(m,n) = y x (mem, M,n+n,N) (2.36) 


l 
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BE ZSTRUCTURES OF TWO-DIMENSIONAL FILTERS 


1. One-dimensional Digital Filter Realization 


One-dimensional digital filters, in general, can 


mensnown in the z-transform domain as follows: 


We 
L 
F 
) Dy 27 
2=0 E 
H(Z,) 5» ———————————— v——c (229 7) 
I M X (zi) | 
-m 
l- } an 2 
m=1 
or, in the time domain 
M L 
AS A) = ) b, x(n-£) (2238) 
m=1 2=0 


The realization implementing equation (2.38) directly, is 

known as the direct form 1 and shown in Figure 2.7. 
Equation (2.37) can be rewritten in a slightly 

different form by introducing a new variable, W(zi), such 


that 


W (zi) Y (zi) 


= ——— 0 — 2s 
H(z) X) WTZ) (2:39) 


The corresponding set of difference equations consists ol: 
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T Uor ord: Aa DI 


T=u 
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w(n) Zar: ) as w(n-m) (2,40) 


m=1 


and 


y (n) by w(n-2) 252313 


ll 
I œm n 


Equations (2.40) and (2.41) are realized in the 
direct form 2 as shown in Figure 2.8 which can be visualized 
mena Cascade realization of two digital filters, realizing 
the denominator and numerator polynomials, respectively. 
Figure 2.8 can be redrawn as shown in Figure 2.9. The 
resulting realization is called direct form 3 or canonic 
form, since it represents a structure with the number of 
delays equal to the order. 

At this point, an additional realization will be 
we roduced, the direct form 4, which combines the charac- 
Perzıstics of the direct forms l and 2 having only two 
entries in each summer, which corresponds directly to 
hardware implementation and by realizing the numerator and 
denominator polynomials separately in cascaded form 
cure 2.10); 

There are several other methods to realize a one- 
dimensional Me order digital filter. For example, the 
cascade form of first and second order sections H(zi), as 


shown in Figure 2.11, where 
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H(2,) = a Hi (21), k 2 —— (2842) 


Each component H, (z,) can be realized in one of the above 
outlined forms. 

Partial fraction expansion methods applied to 
ExU3bron (2.37) lead to a parallel form realization of 


Est and second order sections (Figure 2.12), i.e., 


k 
) = c+ ] mi), k= ES (2.43) 


Other forms of realization include hybrid structures, 

rer, parallel-=cascade forms, the transpose configuration, 
which can be obtained for all previously outlined structures 
by reversing the direction of signal flow and by interchanging 
all branch and summing modes, wave-digital and continued 


fraction expansion filters. 


The direct forms discussed above and the series/parallel 
arrangement of lower order sections are by far the most 
often used ones in computer simulation and in digital 


hardware. 


2. Two-dimensional Bigital Filter Realization 


Two-dimensional quarter plane digital filters, in 


general form, can be shown as 
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Fig. 2.12. Parallel Form 





) ) b Ex Z 2 
V. a ile 2 
4,70 & 470 
H(z, 12.) 0 — m (2.44) 
M M 
alt 2 = zn 
a 2 
1 - ) a Z1 z 
m,-0 m,-20 l 2 
JL 2 


Wie numerator of the two-dimensional transfer function (2.44) 
can be written as a one-variable polynomial in Zor weighted 


by coefficients which are functions in Z3: l.e., 


L L 
2 1 2, 72, 
N(z),2,) = po VN zi )2, (2225) 
2.=0 2.=0 
2 il 
L, 3 
_ 2 
2,=0 
2 


The same reasoning applied to the denominator polynomials 


leads to 


M, M 
EU E 
DC 1 == Lr ^nm, 2 12, 
m,=0 m, =0 
(2.47) 

= 
= 1- Jo (z.))z_ (2.48) 
u M.m J 2 

[ne 
m,=0 
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To obtain the direct form 3 representation, a new function, 


W(z]2,) is introduced: 


ETE ERRE 


HZ } EN Ze 
Jb 2 X (21:25) W(2],2,) 


(2.49) 


where each factor can be written using equation (2.46) as 


L, wi 
2 
W ' = 
2 ,=0 
2 
(200505 
and equation (2.48) as 
2 
zz) + | (Du, ‚m, (21?)2, Y(z,2,) 
m_.=0 
D 
(2:51) 


ime realization of equations (2.50) and (2.51) is shown in 
Figure 2.13 using compact notation, and Figure 2.14 where 
ES Ez. and mm, 21) are implemented for each (L4 $5) 

and (m, m,) by one-dimensional direct form 4 realizations. 

The general form of two-dimensional causal digital 
filter (2.44) can be rewritten by introducing a new variable, 
W(z]2,) Orem the previous function, then, 

W(z42») Y (2,127) 


H(Z1,24) = mm nm 02.52) 
ace X (ZZ) W(2, 72.) 
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usn factor can be written as: 


De > 
“hy 74) 
Y (2,127) = ) 2. 21 2,  MWiz]),2,) 
2,=0 2,=0 +? 
l 2 
a2) 
and 
M. M 
E ns 
W(z],2,) =- X(z),2,) + ) ) Anm. 24 2,  Wiz,),2,) 

LM S 

I 2 

la (2.54) 


The equations (2.53) and (2.54) can be written as a differ- 


ence equation as follows, 


> 
y(m,n) = ) ) P, 2, W(m=2, ,n-2,) (255) 
2.=0 2.=0 
T 2 
and 
n: £5 
emo = xmn) + ) > Emm, w (m-m, , n-m,) (2.56) 
m4 —0 m,=0 
m/m, z 0 
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Without loss of generality, Li is chosen equal to M: for 
EU. The structure realizing equation (2.55) and (2.56) 
is shown in Figure 2.15. It is noted that the number of 
delays in Figure 2.15 equals the order of H(Z,,2,). The 
structure is, therefore, by definition, canonic. 

Inedec £orm 4 (Fig. 2.10) realization of H(2,2,) 
MoesiOwn in Figures 2.16 and 2.17, using compact notation and 
one-dimensional direct form 4 (Fig. 2.10) implementation of 


each N and D , respectively. 
Dita 215 


E STABILITY 

Two major problems in the design of a recursive filter 
are approximation and stability. Since the output values 
are used through feedback by the recursive algorithm, it is 
possible for the output values to become arbitrarily large 
independent of the input. A filter of this kind is said to 
be unstable, which is an undesirable condition. Thus we 
need to know what constraints to put on the recursive filter 
coefficients, so that the filtering operation will be stable. 
This problem is similar to the stability problem for the one- 
dimensional case, except that the added dimension inherently 
increases the complexity of the analysis. The major diffi- 
culty is that the concept of zeros and poles does not hold 
in multi-dimensional systems so that simple algebraic systems 
extrapolation of one-dimensional results is not obvious. 

. Stability Theorens: 
Theorem 1 (Shank's theorem) : A causal recursive filter 


with the z-transform H (2, 127) = (21525 BZ 1120)" where A 
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and B are polynomials in z] and Z>1 is stable if and only 


ENEEhere are no values of z, and z. such that Biz = 0 


i 2 
for all z} > l and z, > 1. (In general we assume 
zB E C NNNM 
I 027009172 DOM M 


that is an expression in increasing powers of 27 and Z3.) 


In other words, the theorem says that if there are any 


1127) 


values (real or complex) of Z and 2, Ioa e B(z,,2,) 
is zero for zz and z simultaneously greater than or equal 
to unity in magnitude, then the filter will be unstable. 

It is much more difficult to determine the stability of 


two-dimensional filters than of one-dimensional filters. 


In the one-dimensional case, it is only necessary to locate 


finite set of roots in the z-plane. In order to mechanize 
the theorem, we must, in general, find the values of z. and 
z^ for which B(Z, 127) = 0. There is no technique for locating 


the zeros of a general two-dimensional polynomial in (21,25). 
Shanks, Tretial and Justice, in their paper [1l] proposed a 
technique to apply theorem 1 to determine the stability of 

a two-dimensional recursive filter, H(z1,25), by which the 
unit disk in the Z} plane must be mapped into the Zo plane 
by solving the implicit two-variable denominator polynomial 
ong H(2,,2,)- A necessary and sufficient condition for 
stability is determined if the map of the unit disk 

EE - (2 


Ik 


d, = (2,,2, che Zo plane. This method requires an 


infinite number of mappings and, therefore, can not be 


112) Peeeecdecs aot Overlap the unit disk 


applied exactly. 
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Huang [20] simplified Shank's stability theorem 


considerably by stating the stability theorem as follows: 


Theorem 2: A causal filter with a z-transform 
Dozen Zk) 
2,2 
Hz) 


where A and B are two-variable polynomials in 21125, is 


Eure if and only if: 


l) the map of the unit circle sd, = (Z, , 12, |=1), 


EScording to B(z,),2,) = 0, lies outside of the unit disk 
d, = (z,.|z,|>1) in the z,-plane, and 

2) no point in the unit disk d, - (21, |z |21) 
maps into the point Z, = 0 by the relation B(Z +25) = 0. 


In order to apply this theorem as a stability test, we 


have to map the unit disk öd, = (z, , 12, |=1) into the z 


l 2 

plane according to B(z,,z,) = 0 or z, = Reise and 

to see whether the contour in the z plane lies inside the 
unit disk d, = (z,,|z,|<1). Also, it is necessary to solve 

8 (24,0) = 0 to see whether there are any roots with magnitude 
Breater than 1. 

Although theorem 2 is much simpler than using Shank's 
Emsginal theorem, the procedure is still infinite in the 
sense that, in principle, we have to map the unit circle 
od, aneo the z,-plane. 


Ansell [21], whose main contribution is to couple the use 


Of a Hermite test for checking stability with a series of 
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Sturm tests checking positivity, has reduced theorem 2. 
Although Ansell's results enable us to test Stable ta 
Ace number of steps, it still is, unfortunately, very 
tedious. 


Make the change of variables: 


- Z 
Py 9 1 + 27 
and, 
7 2 
pem + z 
and let 
DES) 
er ae 
p — — on 
AAA 


where E and F are polynomials in Py and P3- We can restate 
Theorem 2 as follows. 
Theorem 3 (Ansell's theorem): The causal recursive 
filter H(Z152Z») is stable if and only if: 
1) for all real finite w, the complex polynomial 
um "P F(jw,PpJ,) has no zeros in Rep, ? 0, and 
2) the real polynomial in Pi F (p, 71) has no 
zeros in Rep, > 0. 
This theorem is essentially the same as Theorem 2, but 
an advantage of Theorem 3 is that condition 1) can be 


tested using standard techniques of circuit theory. 
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Anderson and Jury [22] modified Huang's method by out- 
lining a procedure which replaces the bilinear transforma- 
Mons by the construction of a Schur-Cohn matrix and checking 
for positivity of a set of self-inversive polynomials. It 
also replaces the Hermite test by a Schur-Cohn matrix test 
and requires a series of Sturm tests or equivalently, a 
system of tests establishing the roots distribution of a 
polynomial. These modifications represent a substantial 
reduction in computations as compared to Huang's method. We 
will discuss stability considerations which are related to 
our design technique in Chapter IV and Chapter V in detail. 

In this chapter, basic properties of two-dimensional 
digital filters have been reviewed and summarized including 
semi-causality to establish a fundamental theory of two- 


dimensional digital filters. 
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III. MATHEMATICAL BACKGROUND FROM ANALYTIC PLANAR GEOMETRY 


In this study, we propose using semi-infinite planes to 
approximate a given frequency response of a two-dimensional 
analog filters by extending one-dimensional semi-infinite 
straight line approximation to two-dimensional cases. For 
example, as we will see, the frequency response (dB vs. log- 


frequency) of the separable canonic factor 


a 


our Ben. un 
npud nu 


can be approximated by a collection of planes. This is 
ENScussed in Chapter IV in detail. 
To prepare for the discussions in the next chapters, 
it is worthwhile to review and develop useful equations 
from analytic planar geometry which are going to be used 
throughout this study with some examples of how to use them. 
The properties of the analytic planar geometry can be 
summarized as follows: 


l. The general equations for a plane are given by 
Ax + By + Cz +D = 0 (O) 
where A, B and C are the direction numbers of the normal 
to the plane. The equation of a plane passing through the 
point (X4 5Y4 721) and perpendicular to the line which has 


direction numbers of A, B and C is given by: 
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A(x-xi) + B(y-y4) v C (z-21) = 0 (372) 


Ene Girection cosines of the normal to the plane 


are given by 


A 


COS Q = (3.3a) 
(A? 4 pg? 4 c^) 7? 
B 
cos § = ——— > ($2535) 
(A^ 4 B^ 4. c^) 7? 
C 
cos y = 
(A* 4 B^ 4 c? 17? 


where a, 8, and Y are the angles that the normal makes with 
the x, y, z axes respectively. 


The relationship between direction cosines is given by 
Cos ^a > aggre + Si = 1 (S226) 


SE a line joining points (X] »Yı +27) and (X,,Y,,2,) 


the equation is given by 


Cos qa o> B = COS Y (3.4) 

n2 1 Y2”7ı e 
4. The equation of a straight line with the direction 
numbers a, b, c and passing through the point (x4 /Y1 721) is 


given by 








sl Zn ei 











5. The angle 9, between two intersecting lines with 


direction numbers a, b, c and a', b', c' respectively is 


given by: 
Cos 0 = Cos a Cos a' + Cos 8 Cos 8' + Cos y Cos y' 
(3.6a) 
ør 
do ar mer ea tbb tec) (3.66) 
ano Aa os 
6. For parallel lines, the relationship between 
direction numbers is given by 
a 2. b . Cc 
at = prt a (3.7) 
ror perpendicular lines: 
aa' +bb' F ccol mF OQ (30793) 


8. The direction numbers of the line of the intersection 
between two given intersecting planes, with direction numbers 


pou C and A', B', C" respectively, is given by 
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a = (3.9a) 
B' E” 
E A 

P (3.9b) 
ES E 
A B 

c = (3. 9C) 
Zt B' 


The angle 8 between these two planes is given by 


Cos 8 = NEU uro er c 


n^ FB eos Au Pal el 
9. The equation for a plane in terms of its slope can 


be calculated as follows: 


Remisceon from Fig. 3.1, the slope, m, of the normal 


line is given by 


m = tan$ = z ay en EN 
A +B? 
m = slope of the plane 
nac rtan(909.—9). - = I = - £ (3.12) 
— an h m . 


12 








normal to olane 


normal to plane 


Fig. 3.1. A General Plane Configuration 
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The equation of the plane is given by one of the 


Following: 
A(x ~ x,) +By + Cz = 0 CORED 
PEN Bly - yy) EZ =O 5. 18b) 
Dr OBL C(z - 23) = 0 3830) 
and 
Ax, = By, = Cz, = -D (3 PSG) 


The intersection of the plane with the z-x plane 


mo Given by 

A(x - xi) FeO T= AO (3.14) 
and the slope, Mr of this line is 
A - $ (3.14a) 


Similarly, the intersection of the plane with the 


Z-y plane is given by 


B(y-yj) *Cz = 0 (3.15) 
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and the slope of this line, My, is 


2 dz B 
E xy - € (3.15a) 


Supstituring (3.15a) and (3.14a) into (3.13) 


yields the following alternate equations for the plane. 


-m(x - xi) Be) ee = O (3.16a) 


| 
o 
W 


-mQXx - m (y - y) + Z 2165) 


| 
o 
C2 


-mX - mjy t (z - zi) 16c) 


The direction cosines of the normal in terms of 


m ram 
E^. o 


Cos a = (3.17a) 
2 2 
my + m^ +] 
Cos B = (3.175) 
2 2 
My +m, +1 
il 
Cos Y = === (3S E76) 
2 2 


Substituting (3.14a) and (3.15a) into 9-11) and 


EZ) yields 


m l (39183) 





and 


2 (3.18b) 


The angle, 6, between two intersecting planes from 


EDU, 13.142) and (3.15a), is given by 


Cos 9 = SES 





2 12 2 


1 +m, +1 


The direction numbers of the line of intersection, 


from (3.9), (3.14a) and (3.15a) are given by 


deve (m, - m, ) EG" (3.204) 
b = (my - my) EE" (22205) 
ec = (m) m, - mom, ) Ce! (3.20c) 


The direction cosines of this intersection line 


are given by 


m5 - m, (3.21a) 


Cos a 





i (3.215) 


Cos ß 
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(120236) 





10. When two planes are added in one dimension to form 


a new plane, the following results. Given two planes 


-mX = My + Z = 2 s 0 
=m! = ' p $ = 
my X may + 2 27 0 
Their sum in the z-direction with z" = z' + zis 
given by 
== = ' We 3 $ = 
(m, +m')x (m, + m, y + 2 (zi +23) 0 


The slopes of the intercepts with the zx and zy planes 


" = - (m, + my) 


3 
te 
| 


" = -(m, + m,') 


3 
N 
| 


The slope of the normal is given by 


m" = [(m, +m1)" + (m, «mj ^j 1? (3722) 


pn 








To aid in comprehension, two examples follow: 

Example (1): Given my (slope in z,x plane), ma 
(slope in x,y plane), and Xy (the point at which x-axis goes 
through the plane). Find the equation of the plane. 


The general plane equation (3.1) is given by 


AX F BY + CZ? tD = 0 
When z= 0 
AY” A 
m = dx _ -B (given) 
3 dy A 
and 
Sap ts (given) 
I A 
When y = 0 
NERA 2 
E C E 
E cA j 
DE Er C (given) 
and 
= . D 
A E 
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If we substitute these into the general equation, the result is 


A om 
Ax Am.Y aoe AX, = 0 
1 
Dr 
X- my ~ er =- X = 0 
3 m 1 
1 
or 
Z = mx- mm,y ~ xm (3823) 


As shown here, we can derive the plane equation easily 
from the given specification. After finding this equation, 


we can find the canonic factor by simply letting 


x = log Wy 
y = log Wo 
Z = 10 log ip}? 


Alternately, we can go from the canonic factor to equations. 


Maus equation (3.23) can be written as 


i : 2 
10 log |T(3u4 7305) | m, log u, - mm, log u, - 10 log K 


(3.24) 


where 


10 log K = m4 X4 











This result assumes that 


2 x 
IT (394 7305) | = E (524a) 
q 
Ku, 





Making the logarithm of (3.24a) yields 


2 
MENOS |T| = 10 p log o, - 10 q log w - 10 log K 


(3:525) 
Byeequating (3.25) and (3.24), we get 


10 p = my (3.2023) 


10 q = m (3.26b) 


Beecoliows from (3.26a) and (3.26b) that 
m, = Qq/P - 
mas the transfer function of (3.24a) has a logarithmic 
characteristic as sketched in Fig. 3.2. 


The previous result can be extended to the following 


transfer function 


(1 + ju,)?/? 
uu 729 ) — EMI IEEE 72 (955209 
l E xt A ~ w = 


Besubetituting (3.264) and (3.26c) into (3.27), then 
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Eg. 


3. 


2 


Dres logarıctımrezcharaetreristie of 





m p 
* . 2 PM i 
|T (jw JU ) | = 
1 2 q 
Kw 
2 
Qi kee; au — m, log w - mm, log 6 
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2 


_ CONES 





T(Jw, ,Jw),) 


when w 


T(joi 73u9) 
and 

20= log TN 
When w << l and D Jie 

T(ju4 Jw) 
and 

209g mi 
When DS Mand Dec l, 

T(ju4 Ju) 
and 

to | 


| m4 /20 
(1 + Jo) 
m 4210 

«1? (1 E ju, 2 


(3. 28) 


m l and Wo << I region I 


x 17? 


- 10 log K 


region II 


Tm 


me 720 


g1/ 2 


M 


= IL) ee V ane doct 


2 2 


region III 
| m,/20 
(jo) 


„1/2 


O A ch kE my log 4 
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when Wy >? once XT region IV 


2 
| m,/20 
Jw 
T(jw 2555 ) = 
1 2 VA m>,/20 
K (Jus) 
and 
20 log |T| = - 20 logK +m,log u, - m, log u,. 


The resulting planar approximation for the logarithmic 
transfer function is shown in Figure 3.3 
Example 2: 
In this example, we derive planar equations from the 
given transfer function. 


Given a separable two-dimensional transfer function: 


IM" s K 
T(jui»205) =  —__———— (3.29) 


(1 + ju, 11) (1 + ju, c2) 


The magnitude square of the transfer function (3.29) is 


given by 
2 
2 K 
T(Jw¿rJu )| A  _ __.—__u———_=> (3:30) 
| Sa (La 0 (140,1, Y 
Make the logarithm of both sides of (3.30), 
ae loa ogla og 1 +20. 109 R 
10 log |T| = -10p log(1l+0, "17 q log > 79 
(0379 
when 0414 << l and Wats << Jl Ecnoeceguatuonei3. 3l becomes 
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2=2ole ITI 


neben, 















q = logo 
Ys log v^ 
> 2 log wy 
Luc 3.3. The logarithmic characteristic of 
DE 
Tues D c 
il 2 ki/?( F ju, 12/20 
-10 log K + m, log Wy 7 mM, log Q^ 


20 T 
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2 
ON En = 20 Log re 


where C is some constant. 


When Q4T4, >> l and Wats >> L; 


2 
Ur loge T| = -20p log wt} 7 20q log At, 
+ 20 log K 
or 
2 
10 o |T] = -20p log w, - 20q log wz 7 20p log T} 
> 20 10g tT, + 20 log K 
Let 
2 
Zee ee LOG hr | 
Ks Log 03 
y = log "D 


and we get the planar equation, 


2 —20px — 20qgy =- 20p Tog EE 20q log To + 20 log K 


or 
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um 
(-20p)x * (-20g)y - 2 - 20 log ( O 
or 
mx tmy-z - Ky 
where 
m. = SlOpeseu tier linem(x<, 2) plane 
= =-20 p db/decade. 
m, = slope of the line (y,z) plane 
= -20q db/decade. 
In a similar manner the regions (0111 REDE 
Wats >> and (wıT] =>]: Wats << 1) can be investigated. 


The planar approximation for the logarithmic transfer function 
e sketched in Fig. 4.1 where q = p= l. 

In this chapter we have developed some useful planar 
geometry formulas and shown how they can be used to determine 
a planar approximation for the two-dimensional logarithmic 
transfer function for separable cases. Some non-separable 
factors are discussed in the next chapter in which experimental 


verification of this approximation technique is also presented. 
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IV. CANONIC FACTORS IN CONTINUOUS FREQUENCY DOMAIN 


In Chapter III, we discussed some equations of planar geo- 
metry and showed how to use them for approximating separable two- 
dimensional transfer functions. In order to improve the scope 
of the proposed design technique, we now consider the behavior 
of other separable and non-separable factors insofar as 
their planar geometry approximation is concerned. This will 
provide insight into the frequency response characteristics 
of two-dimensional systems and indicate what can be expected 
by cascading these canonic factors, knowing their individual 
Benavior. The logarithmic characteristic of cascaded factor 
HcwSsumoty the sum of their individual logarithmic character- 
istics. This approach parallels the semi-infinite line 
approximation (dB vs. logo») technique which has been used 


SO successfully in one-dimensional design. 


A.  SEPARABLE FACTORS 


The following several separable canonic factors are now 


considered. 
(1) Be ae le ie) (4.a) 
(2) e Dei) T Da T i (4.b) 
j Did 
(3) pU (4. c) 


where p, q and o are integers. 
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(1) First, as an example of the technique, consider 
the following analog separable transfer function, the first 


separable canonic factor with p = q = -1; 


I 


ME ee 4.1 
p (Fjot) (lF jut) m 


The logarithm of magnitude-square of this function 


is given by 


2 2 2052 
10 log |H(w,+w,)]% = “10 log(1+w,%1,%) - 10 log(1+w,t,") 


(4.2) 


Borsus define the following for simplicity of notation. 


2 
z = 10 log [|H (94 702) | = 20 log |H(w] +09) | (dB) 
(4.3a) 
x = 20 log Wy (4.35) 
y = 20 log (^ (4. 3c) 


Consider the regions in the (w rwa) plane (as seen 
in Figure 4.2) as follows: 
Region I, when wT] < l and WT, < je; 
then z = 0 db (4.4a) 


Region II, when wT] > l and WT) < d 


z = -20 log (a4 t1) db (4.4b) 
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Region III, when WT} < l and Wats > T 


VERE log (w,T,) db (4.4c) 
and 


Region IV, when WT] > l and WT, > D 


z = -20 log (a4 14) EDO (a4 14) (4.4d) 


The above equations (4.4a), (4.4b), (4.4c) and (4.44) 


can be expressed as 


Region I, Zo) (45252) 
Region II, Z = =x = 20 log Ti (4.5b) 
Region III, z = -y - 20 log T, (4456) 
Region IV, IZ OO E T9 7 20 log T1 
(4.5d) 


These equations describe four intersecting planes in 
X,Y,Z space as shown in Figure 4.1. 

Plane I is the horizontal plane at zero db. Plane II 
and Plane III have slopes of -20 db/decade passing through 
the corner lines as indicated in Figure 4.1. These lines 


are analogous to the break point frequencies of one-dimensional 
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Zu 20log Hlw,.w, ) 


Wija WFO e X = log Wy 





slope —-—2048 





Slope = = 2048 


Na Ll. Planar Approximation for 


MR NE. .— 0 en 
ame (1+jw] 17) (1+jw,T,) 
20 log |H(w. ‚w.)| = -10 log (1+w ^ a 
3 5» De 
2 2 
-10 log (Itw, > ) 
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X =log W j 
W = Yn 
by 
Gigs 4.2. The regions of 
H(jw jus) : 
Jt 2 (1*j0414) (1530515) 
20 log |H(u,,9,)] = -10 log (l*u, ^c, ^) 
3 Lee i € 
DN 
0 sep Usa sen) 


in (log 24 ,log Wo) plane 





log-modulus plots. Several properties of this plane are 
summarized. Pertinent equations from analytic geometry are 
given in Chapter III. 

The actual frequency response of this separable 
panonic factor with 1 = 15°51 is shown in Figure 4.3a 
(dB vs. log Qr fog Wo) which fits the planar approximation. 
In Fig. 4.3b, the contour plot of the frequency response is 
given, which verifies the regions in log we log Wo plane 
Shown in Fig. 4.2. The logarithmic plots versus linear 
Maeauency are given in Fig. 4.3c and Fig. 4.3d. It is 
apparent that the approximation regions are well defined 
when logarithmic frequency is used. 

The frequency response of this separable canonic 
factor with E gos tT. = 2 iS given in Figures 4.4. These 
results again confirm the logarithmic approximation. 


(2) Next consider the canonic form 


y 


H(jo4 73023) e (1 ju512) (1 30414) 


where we have two separable factors of "PE The frequency 
response of this factor with LOSS l, T3 = 2 and D 4 is 
shown in Figures 4.5. As seen from this figure, we have 
-40 dB/decade slope in the direction of Wor as expected 

from the one-dimensional case. The regions are shown in 


the contour plot with dotted line (Fig. 4.5c). 
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In a similar manner we can analyze these regions as 


follows: The logarithm of magnitude square of this function 


is given by 


2 


2 
10 log |H(w,,w,)]% = - 10 log(l+u r 1) 


ace 10 Tears 


2 2, 
Wo Tz )- 


- 10 log(1+ 
Consider the following regions in (01 504) plane (as 
Seen in Fig. 4.5a) as follows: 
Region I, when q Ty << l and Wats < Wot << i 
Bene - 0 dB. 


Region II, when "pm d l and Ut € U4T4 << I 
zZz = = 20 log CTESE dB 


Region II, when sy ee Isanda T cT 


ao e 3 


-20 log (w ) =» 20 log (w dB 


N 
le 


272 273) 
Region IV, when q Ty >> e WoT3 > WoTo >? 1 


T 


EE COE 


N 
la 


LIE LA Lec E 


en 2.09 
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Region V, when wT] << E QAt4 >> Es WaT3 >>Woto >> 1 


2 = -20 log (w573) = 20 log (w,t,) 


Region VI, when wT] >> d WoT3 > Too >> l 


Za == 20) LOS = 20) Log sa acd) logo 


l't 2:05 aa 


where 
Z = 20 log [Hwy 119) | 
(3) Consider the separable canonic factor now, 


H(u,u) - [1+ (Ju 11" 


We can show that the following regions occur by 


punc the notation of (493a) sand 4.3b) for simplicity: 
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Region I, when Qi ? 1/14, 


Za = px + 200p fog Ti dB 


Region II, when Q4 < 1/14 


These separable canonic factors are identical to 
the canonic factors of a one-dimensional transfer function. 
They are just an extension to the two-dimensional case. 


Therefore, they need not be discussed further. 


B. NON-SEPARABLE FACTORS 


The following non-separable canonic factors are now 


Poneidered. 
(1) Berger Que) ae (4.4) 
(2) p (legten eese e DEL (4.e) 


where p, q and a are integers. 


(1) First, consider the non-separable canonic 


Be@tersor the form: 


H (a4 (0) = [1 + (joy 11 )P (3319117 
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Region I, when 


ZEN SAEZ OG Ta 720 aloe T dB 


The cornerline between the two regions is given by the 


equation 


Or 
pu scm ODE log Es) zug Tog PE 0 
The cornerline is sketched in Figure 4.6. 
The slope of the plane in region II is given by 


-d/o* a a . Direction numbers of the normal are given 


Bu cd, ag and -l. 
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Depending upon the sign of p and q and the magnitude of 
T1 and T Several different orientations of the cornerline 


are possible. 


Example (1) 


Discuss this canonic form in the special case with 
a= q= l, a = +1, T; = 0.5, and > = l. The transfer 
function gets the following form: 
+1 

H(w, , 07) = l F (.5507) (Jw,)] 

Take the logarithm of the magnitude square. 
2 +2 

|H (w rw) | = |l - 0.5040» | 


A) er = 504 0» | 
UpSregrlon I; Wit, << 1 and Wats << T 
z = 0 dB. 
Eneregrlon I1; wT Sei ana Wats >> Jl 
z=x +y + 20 log (0.3) 


where 


x = 20 log uy, y = 20 log u, and z = 20 log |H(o4 502) | 


The actual frequency response of this factor is shown in 


Fig. 4.7 including a contour plot. 
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The contour plot of the frequency response of 
H(ju4 , 02) m (301) (30.505) 
z = 20 log|H| = 20 log (1 - 0.5040») 


in bod Q4 , log Q^) plane 








When a = -1 we have a stability problem in this example. 
The denominator of this transfer function becomes zero 


when 


O 


I 
e 


12 


Dr 


Wy Wo = 2 


This is the singularity of the transfer function in (04 107) 
plane. It can be called a singular line of the system. It 
is a curve in the (07,07) plane to be compared with a pole 


in a one-dimensional transfer function. 


Example #2: 


In this example, we discuss the same canonic form 


Ene different coefficients, i.e., 


Then the transfer gets the form of 


25 
u 1, +1 
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Taking the logarithm of its magnitude squared 
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Eu su) |” = [1 + —4/? 
n 
2 
901 
2 — ZO Ol | 
(304 
In Region I when Sc 
22 
Z = 0 dB 
307 
In Region II, when " >> I 
2 
Z = 20x - 20y + 20 log (.5) 
where x = log u), y = log u, and z = 20 log |H(w,,w,)|. The 


boundary between these regions is given by the equation 


We also have a stability problem for this canonic factor when 
u, = (mais is also a singularity in (w rwg) plane in the 
form of a straight line. This type of factor is unstable. 


(2) Now, consider the non-separable canonic factor 


u p.a 
H(w, 107) = [l + (Jw, + juzgua) ] 


We have two regions again, as follows: 


Region I, when wT] + Wats <$ a 


AR > 








Region II, when WT + Wot > 1l 


z = ap log (0111 * 051.) 


The separation line between the two regions is 


given by the equation: 
Q1Ti + wT) = I 
Or 
log (a4 1 t 9412) = 0. (4m6) 


The locus of points of (4.6) is shown in Figure 4.8. 
Borsreglon Il, we have tozthink Of two Conditions, 
as follows: 


a) When wT} > Q4 T^ 
zr og Ta 
b) When wT < WoT5 


ZEE oy. r ap20 10g T4 


Thus, region II is split into two areas with the separation 


line given by 


E Bo 


TL6 








Ww. sido wis x=logw 
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a c ec X 


Y-log vw, 


HEEL 428. Locus points of log (wT, E Wats) = 0 


IET 








Frequency. response of this non-separable transfer function 


TE ovn in Fig. 4.9 with p = 1, a = -1, T1 = T, = l. 


ee CONCLUSIONS 

In this chapter, we discussed the properties of the 
separable and non-separable two-dimensional properties of 
canonic factors inthe analog domain (07,0, plane). We can 
conclude that these canonic factors can be approximated by 
planes in the log magnitude versus log Q4 and log w, Space. 
As a design technique we can obtain any given frequency 
specification by proper combinations of planar approximation 
Of canonic factors. It is obvious that only certain special 
cases and symmetries can be obtained using the factors dis- 
cussed in this thesis. For example, a wedge-shaped character- 
istic in frequency domain could not be synthesized by means 
of the factors we have considered. 

The other result of this study is the stability problem 
mach arises in certain transfer functions, that is, the 
singularity curve or line. For a separable case, we don't 
have a stability problem; but in non-separable cases we may 


have stability problems. 


12158 











— 
ee 


V 


a 
\\ 


\ 
IN 
| 





Y 
\ 


7 
\ 
Ut 
i 


- EC 
A e 
SNS 
AS Ss 








a 


— 
—— 


-enne 


\ 
As 
Y NS Ek EO 
NS N s Lom 
\ SS x —— 
N 


—À 


N 
N 





RU 












EEE 
©, 20 de o 
O QA 

Soo 


C05 
5 
3 % CAI 
0% s 








pp 2 
EEE 
(7 LTE 
A RE Q5 
PLL LIS AS 

EEE 
LOO LACH G 
dy} Mee I e 
f (/ Y N n», yy UP on = 
HH N 
// 















| | 
L 




















| 


I 












M», 
2 





Ms Y YY 

// VY, MN N o os = 
TU M —  — ——  — 
7 — —— 


HE 
oe 
O y = 


m =E 
Li HUI HY 1 a 
| IM z= 





i 





HM 
5 





4.9a. 


Erg. 


l 
Jr JU, zi ju» 


The frequency response of 


1 


) 


rju» 


i 


H(jw 


— 


tus) | 


20 log|H| = 10 log |1+(w, 


Z = 


(dB vs. log u,,log w,) 


129 





l mboy 


OQ! 
| 
} 


"0I 


logw , 


Ac CHIESE SENSE NC. 
f=sCHLE=1 SOUR dame eS INCH. 


Fig. 4.9b. The contour plot of the frequency response of 


GET z 1 
Duo UNE I el * ju; 


z = 20 log|H| = 10 log |1 + (wtw) f| 


gm LOG Q4 ; EOG Q^) plane 








Q 


e 
(2 


Y 
$ 
X 


(à 


Y 
O 


(à 


» 
$ 


o 

Y 
Q 
x 


Y 


e 
x 


x 


(à 


Y 
(à 


Y 
Q 


o 


0 
0 


X 


Y 
© 


Y 
o 


Y 
N 
è 


N 


Y 
Q 


A 

% 
K 
X 


$ 
X 


Q 


Y 
o 


Y 
^ 


Y 
N 


" 


N 


O 


Y 
^ 


A 


Y 
O 


Y 
2 


Y 


Y 

Y 
Y 
x 


| 
\ 


Y 
NN 


Y 
. 


Y 
o 


Y 
Y 
N 


X 


vi 


Y 
Y 
N 


Y 

Y 
Y 
N 


(à 


Y 
x 


0 
A 


e 


N 


x 


Q 


Y 
A 


Y 


Y 
Q 


m 
X 
X 


Y 


Y 
0 


Y 
N 


Y 


X 
A 
Y 

N 


O 
A 


Y 
N 
" 


» 


7 
X 


Y 
X 
X 


Y 


O 


^ 


ae 
SYS 
N 
3 


Y 


Y 


Y 
Y 


Y 
N 


E 
A 
i 


E 


x 


o 
Y 
XX 
NN 
© 


V 
Q 


Y 
X 


Y 
Y 


Y 


v 
9 


Y 


Q 
O 


i 


V 


Y 
S 


Y 
Q 


x 
X 


5 


Y 
Q 


S 


X 


X 


Y 


Y 
Y 


Y 
Y 


v 


Y 
O 


Y 
Y 
> 


Y 
Q 


Y 
Y 
N 


ok 


à 
S 
Q 


D 
Y 


O 
O 
x 


2 


Y 


















2 OO OR 
OLD (Se $ SY X OY $ LOY OX AX AS OX SQ Y QU 
SEES E SS RS RE 
OS SOS LS PK RR OG LEEREN 
RS x SS OCR RRO LER EEK ET AYN 
ER ^ SS QO AX ROR RRO RCT RER RN: 
SX > RR SS ¿2 S EQ ene AQ ER OR OR AID RI 
S SÍ ISS OR SCO Hy I IG es RRR KORE 
SR RS SS SIA RR 
> e O SS OX CS oY © RR A S COO ORO S Qro CODO 
ES RS SS RR Xe ERLERNEN 
Se SR RT ERIK Se x REKEN SO ER EN, OR 
SS SS ES RAS, A ER OR 
694029494500 0 054549, O OR A ROSA RE 
SS SS o 9e oS OS ES ER X) XO RES LTE, LER QU 
IR S EEE SSO MRR ORI N 
QS, ISR oS eS ES QS O X Y e AS ER OUS ET O " A SL ^ ^^ LEN 
SL RER EEE IN I I IK COSA KAR KER MAK RR) 
ERROR 9,9, 4 0 4 05 XC S II RO t RO SEI 
SEES OR O SOI SR NA 
SÓ RRS KKK OKC IR IH III eA KICK 
SISSI O O IO O RU hey 
CERRO O DOSIS OSI AAA 
A SEESI II OS ROI OS 
SIR LR ORY LLO RR OS A RO OO 
wo IR RR LEERE EEE AAA 
OOOO EL SS Sg ODO ROSSO 
RTL IL AOS AS 960 X OX 8 060606094 Y 0€ 6 OLX D Y ete SS NOS NL A 
E RER SESSI RIINI IG OO KO 
res 24954460: EEE ET Ran O A IO O 
RL AAA, 
ae POSE OXON 0% O OO ROCOSOS RR SAI 
RR OR OO II 
¿SOS O RR LS RI SA SAN 
SOS DOC CÓ SRN RRR KEK RRR ORRIN AA II) RR 
OS O AA DRA RENN RR 
€909900006090090600906000900009000060000000100000009 0.0 COO DOGO 
SEO LORI SAK K RX MY OR RANN A RAKAA ONAN MA 
O 04000050 D 299 $ 
SLRS DIO 0 0,0 0,9 0 99 606 6 05 6 QOO OR AA AS 
III KK OO) QOOOOOOOO 0f 
C OX OX X Q0 060€ 0600 e PLLC KIKI IH IK 0 n 
0050... .0 0 0 0 0 0 0:0 0.9 0,0 /0,0,6 40,0, 9, 5,5, et S tite SI 
SSR SSO XY) RR ee eR OCT 
COS SRL KK RKO SN III 
SER CSS AO SO INSI A NS AIN AI 
060,960 69606060 060 85 RR LO OO SOSA AAA RR 
SERRA 2525509 OIR ROO O 09007996200 SNANAR 
IR O CORRAN 9 9 9 0 990 90 0 069690595080, (26959,9289 90) 90 90 OO NK RY 
ee 0S0 S00 05.5.00 0 000€ 060, 99050596 90505 9 9 X EEE EVENT 
s SOS 252505 70:0 0.0, 0,0 ICICI III ODO OSOS OSI III 
SSSR) SO BRK IKK III KK ARA AO 
Sd, IO O A 
SLS LSLS OR RR SR IO IO OA 
SE RR ERROR RIO OOO OO OSOS OOOO SR IO A) 
SERRA A AA 
C X2 SR X RR OO SNS IA RO A ANY 
LIS ZSERIES IIRL RR RS GTI AN QUE 
C2525 9.0, SLI DO LLO os RS e. QS QS 5 OX XN XS SS AO OA SA 
Qo 0060050 2009000 X ER BKK “ e OR ER o 
IY 


+ 
è "n QOOOOOQQQ OSSA y 
BK CKO KG 
RR AAN 9) 





2 
(linear frequency scale) 


1: 


e 
2 L + Jv + Ju, 
Q4 02) 


Jwy ,J0 
121 


20 log|H| = 10 log |1 + (w,tw,) 


The frequency response of 
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V. TWO-DIMENSIONAL DIGITAL FILTER DESIGN TECHNIQUE 


A. INTRODUCTION 

The design of two and multi-dimensional recursive digital 
filters is difficult due to the absence of the Fundamental 
Theorem of Algebra in two or more dimensions [15]. That is, 
polynomials in two variables may not, in general, be factored 
into lower order polynomials. For this reason, several design 
procedures for two-dimensional recursive filters have been 
proposed in the literature. We then consider the use of the 
double bilinear z-transfer with the semi-infinite plane 
approximation in the previous chapters. These are briefly 
discussed in the following. 

Shank et al suggested three design techniques in their 
paper [11]. 

l. The use of a series of stability-preserving transforma- 
tions to produce a stable two-dimensional recursive discrete 
filters from a stable one-dimensional continuous filter. 

2. A space domain design technique in which the unit 
sample response of the filter being designed is made to 
approximate a given desired unit sample response. 

3. A stabilization technique in which the denominator 
polynomial is replaced by its double planar least-squares 
inverse. Planar least-squares can be defined as follows: 

We have some given array C. We would like to find 
an array P such that C convolved with P approximates the 


unit pulse array U. That is, 


C*P zx U 


1022 





where the symbol * denotes two-dimensional eonvelutson. In 
general it will not be possible to make C * p exactly equal 
to U. In actuality, C * P will be equal to some other 
array, G. If we choose P such that the sum of the Squares 
of the elements of U-G is minimized, P is called a planar 
least square inverse of C. Planar least square inverse of 
P is called double planar least square inverse [11]. This 
new polynomial is conjuctured to have no zeros on 


—2 
u = 1{2,,2,): E >,|z,| > 1} and to have a frequency 


I 


response which closely approximates that of the Original 
pronsrer function. 

The first technique is known as a filter design with 
rotating axes and these filters are called rotated filters. 
This design procedure consists in considering a transfer 
function in the analog frequency domain, rotating the axes, 
and applying the bilinear transformation to obtain a two- 
dimensional recursive filter. In general, the starting 
point is a separable filter and stability of the resulting 
filter is not guaranteed even if the prototype filter is 
Stable. This technique was used by Costa and Venetsanopoulus 
[23] to design low-pass filters having circular symmetry, 
and also they proved that rotated filters are stable if the 
angle of rotation is between 270° and 360°. A technique 
similar in character but rotated axes in digital frequency 
domain ((z1 525) domain) has been given by Chang and Aggarwal 
[24]. Hirona and Aggarwal [25] have accomplished the rota- 


tion through the introduction of rational powers of Zz} and 
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Z4, and the interpolation of input and output sequences. 

They have proved that, in this technique, the resulting 
filter is stable if the prototype filter is stable. Further, 
the distortion, a consequence of the bilinear transformation, 
is absent in the resulting filter when compared with Shank's 
Señal. 's filter. 

The authors of [11] were unable to prove that the new 
denominator polynomial produced by the third technique, 
which is an extension of the one-dimensional case, was 
stable in the two-dimensional case. Nonetheless, several 
hundred filters have been designed using this technique 
without a counter example appearing [18]; furthermore, Jury, 
Kolavennu and Anderson [26] have shown it to be valid for 
certain low order cases. However, Genin and Kamp [27] have 
recently shown that counter-examples do, in fact, exist. 

Spectral factorization is another technique which can be 
used to produce a stable filter from an unstable one. By 
Spectral factorization we mean a complex map that carries a 
Stable rational transfer function into another stable transfer 
function exhibiting a different frequency response, at the 
same time maintaining some desirable characteristic [29]. 
Pendergrass, Mitra and Jury [28] used stability-preserving 
spectral transformations to modify the frequency response of 
quarter-plane filters. 

Mitra and Chakrabarti [29] extended this work. Read 
and Treitel [30] suggested a two-dimensional Hilbert trans- 


EO perform the spectral factorization, while Pister 
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[31] and Ekstrom and Woods [32] made use of the two- 
dimensional cepstrum.  Ekstrom and Woods also suggested ways 
that the factors could be truncated and smoothed to preserve 
stability. They also showed that the usual quarter-plane 
filter is not the most general two-dimensional recursive 
filter. It is instead possible to have a recursive filter 
which in the limit has a unit sample response with support 
on an asymmetric half-plane, and also they developed 
stability conditions for this new class of filters. 
B. TWO-DIMENSIONAL DIGITAL FILTER DESIGN PROCEDURE BY 
DOUBLE Z-TRANSFORM 
Double z-transform 

In the one-dimensional recursive digital filter 
design, the bilinear z-transform is often applied to an 
Analog filter transfer function. We proposed that this one- 
dimensional design approach can be extended to the two- 
dimensional case. Before going into the design procedure, 
two-dimensional recursive filter design via the double 
bilinear z-transform, it is better to prove that we can 
use double z-transform for two-dimensional Systens. 

Each of the canonic factors, proposed in Chapter 
IV, either separable or non-separable, can be converted 
into discrete equations by the double z-transiorm, aS 


follows: 


(5.la) 


125 








i -1 
2 *2 (5.1b) 
S — een fee: : 
E B dU t 


2 


The single z-transform is equivalent to trapezoidal 


integration which can be shown as follows: 


X(s) = S Y(s) (522) 
x(t) = a 
OF 
Ti 
DUE = EZ clie 
nel 


) (534) 


I|: 
| 
s 
+ 
A 


where T is the sampling period. 


We write (5.4) in the z transform domain as 


ü-zlvcz) - Züasz x) 
Hence 
2,1 - a (en 
X(z) = (=) Y (2) : 
IT 12 if 
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eemzarıng (5.5) with (5.2), it is seen that 


21 
2 l1 - 
s = 3 =) (5.6) 
JEDE 
Equation (5.6) is known as a trapezoidal bilinear z-transform. 


We can drive the double z-transform in a similar 


way, as follows: 


a? 
W(x,y) = dxdy 2(x,y) (S274) 
er 
W(s, 57) = 5,5, Z (S4 S3) (5 I5) 
Then 
Í fa? Tl — EID DNE EC (5.8) 


by a discrete approximation, the left-hand side of (5.8) 


becomes (see Fig. 5.1) 


ll? 


f fa? z(x,y) P E 1)] 
[2 (x, 7 Y4) 5m EY 


-1 -] -] -1 
(1-724 - Z5 +2] 2, )2(2,12>) (5.9) 


127 








and the right-hand side of (5.8) is approximated as follows: 


{[{W(x,y)dxdy = SL ,y ) DU v 


2E ecc) p ESQ 


_ AxAy Sa =I holt 
mu = 2] - Z> + 2], 2, JW(z,,2,) 
5.10) 
By equation (5.10) and (5.9) we have 
7-7.) 1 - 271 - Eu ~ 2. 12-} 
p 2 A 1 2 1 2 4 al 
Z2 5) -] =] -] -1 ` AxAy : 
2272 1 + zi + Z > -— Z1 2, 
GI 
1 - a 1 - se 
ne, 1, 20. — 2.) (5.12) 
AX -l ^ Ay -1 i 
l + Z} l + 2, 


BEconparing (5.Ia) and (5.25) it is seen that 


-1 
l- 2 
S = E 1 ) 
n AX -I 
l + zZ 


[A 


S = HA 
2 Ay l + Zo 
O BD. 


Thus, this is called a double z-transform. The 


double z-transform seems to be a reasonable substitution 
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y (1) x (t) 





(a) One-dimensional case 





(b) Two-dimensional case 


O DL Two-dimensional extension of the 
Bilinear Transformation 


152.9 








for each canonic factor. However, each factor must be 
checked for stability, which is discussed in detail later in 
emis chapter. If each transformed factor is BIBO stable the 
overall transfer function will be BIBO stable. But sometimes the 
double z-transform may in certain cases lead to unstable 
solutions [33]. Finally, it is noted that the double z-transform 
results in a quarter-plane symmetry which is a severe limitation. 
2. Design Procedure 

The two-dimensional recursive digital filter design 
investigated in this thesis is the double bilinear z-transform. 
Proceed as follows: 

l. Determine a two-dimensional analog transfer 
iüamsctron, 


H (s1 $5) = B(s,,5,) CS EM 


where A(s,ıS,) and B (S4753) are mutually prime two-dimensional 


polynomials and B(s,r5)) satisfies 


B (S4753) z 0 
(37139 


A 
V(s,:59) e T = ((sy°S2): Relsy) > 0,Rels,) 2 0) 
in such a way that, the frequency response of G(s} S3) 


meets the given specifications (for example as discussed 


in Chapters III and IV). 
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2. Apply the double bilinear z-transform 


E -1] 
D CI 
S, = £~ iia) 
l + z 
1 
ES a 
Ms) (5.14b) 
l +z 
2 


where A is the sampling period, which is assumed the same in both 
directions. This yields the two-dimensional digital transfer 


lunction 


C (21525) 
D(2,,2,) 


A 


HD CD) - Tz) E) 


where C (21525) and D (21,25) are mutually prime two-dimensional 
polynomials in Z1 and Z5. (5.15) will have the desired 
frequency response with proper sampling time without 


aliasing effect. Since 


1 - Fg 
San (—) 
l + 27+ 

3 


Ber cher elosed unit disc onto the closed right halt 
plane, one would expect (5.15) to satisfy the following 


sufficient stability conditions [11], [20]. 


D (24725) x 0 


qos 





-2 A 
Virza) c u^ $ ((2,,2,): EA AN (5.16) 


Unfortunately, however, the conditions (5413) end 055 15) 
as we will show later in this chapter, are not always 


satisfied. 


In summary, we have two conditions to be satisfied 


Ber stability. 


i: B(S),S5) z 0 
Des) er = {(s],S,): Rs], 2 O0,R.S, » 0) (5.17) 
Ele D(z2,,2,) X 0 
Vc A) E y? MA RN A A E 1 
1 = 2 Jo ro qu EE 210 


(Salon 


Meme Y stands BOG eros each”. amd u^ is the closed unit bidisc. 


C. DIGITAL CANONIC FACTORS AND STABILITY CONSIDERATIONS 

D xmcsductrion 

Certain properties of two variable polynomials and 

rational functions which will be needed later will now 
be discussed. It is well known that a two variable polynomial 
is not, in general, factorable into first-order polynomials; 
but, a two variable polynomial can be factored into irreduci- 
ble factors which are themselves two variable polynomials 


but which cannot be further factored. Of course a given 
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polynomial may itself be irreducible. These irreducible 
polynomials are unique up to multiplicative constants [33]. 
Two polynomials which have no irreducible factors in common 
are said to be mutually prime. Consider two variable 


rational functions, in general, 


P (21725) 
tal T aapa 


where P (2,127) and Q(2,127) are mutually prime. Using 
terminology of [34], a point (2, 127) making Q(z7 125) = 0 
But P(z,12>,) #0 is called a pole or a nonessential singu- 
larity of the first kind (such a point is analogous to a 
pole in the one variable case). A point (23,25) making 
P(z],25 = Q(2],25) = 0 is called a nonessential singu- 
larity of the second kind (such a point has no analogue in 
the one variable case). Clearly, if (27,25) is a pole, 
G (21,25) = ©», but if (27,25) is a nonessential singularity 
of the second kind, then G(z],25) is undefined. 
Zee parable Factors 
The most common canonic separable factor in Chapter 
IV is 
Bu un a en rn... (5.194) 
Ir 772 


(1 + ju4 11) (1 + JW5T5 


assuming Ure TO Pou Stitt timc 


P3 





We get 


1 


Ev S " (1* si) (1* s.) 


Substituting equations (5.1a and 5.1b) into (Sie oD), 
performs the double z-transformation. Thus, 
8127) (1 z;l 


Z) 2 
I 2 (a +277) (a z51) 





2 


where 





A? 
p elo 
2 
a = 242 
ao 
A = sampling period 


5.190) 


520 


This separable canonic form obviously satisfies the condition 


of (5.17). We should be careful about the second condition 


ER (5-18). 
When Z, = 25 = -], we have a nonessential singularity 
of the second kind if a = 1, but a = 1 never occurs. 


Therefore, there is no nonessential Singularity of the second 


kind. Thus, this separable case is stable by satisfying 
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two essential necessities (5.17) and (5.18). The frequency 
response of this separable case is plotted in Pug. 5.2, With 
sampling time A = 0.08 which fits that of the analog frequency 
response. 


The double z-transform for equation (5.19a) is 


2 s lı 1+ 2-1 
T(z_,z_) E EEE e ee 
D 2 (A-21,) (A-21.) A+27 A+27 
1 2 (— 27% ( A 
A-7214 l A> 21, 2 


and this general case for the two-dimensional separable 
factor is stable by satisfying the stability necessities 


Sao. 17) and (5.18) 


The second analog canonic factor in Chapter IV is 


| ! = gl 
H(ju],Ju,) = (l*jo41:4) (1*jo515) (1* ju514) 


When we perform the double bilinear transform for 


this, we get 


A? (1 zi JO 
EZ ,;„2.) = 
072 a A+27 A+27T 
(A-27,) (0-21,) (0-27,) I + 2) 1 +2) 1 2j] 
St 2 neo M y A-2T 2 A-21 2 

1 2 3 

(5.22) 

- 24 = 4 


The frequency response of (5.22) for ae > l, T3 = 


is shown in Fig. 5.3 which almost fits the analog case. 
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Fig. 5.2b. The frequency response of 


i 8(1+2,7%) (142, T) 


(a2, 1) (a2, 7) 





T (zi; 25) 


8 = —— , a = $$ , a = 0.08 
(A-2) 





07107) (linear frequency scale) 
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Fig. 5.2c. The contour plot of the frequency response of 


e a °°) 


T(z),2,) = _— |, where 


(a*z, 7) (a+z, ) 


2 
8 ——, g - B2 , a = 0.08 
(^-2) 


(linear frequency scale) 
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We can summarize the mxn analog separable factor as 


ee a 
1 


H(jw ¿Ju ) = [ 
m a jar j 


-] (5.23) 


loi o 


1=1 


After performing double z-transform, (5.23) becomes 


M A (1277) N A" (1 27) 
T (2,125) = es 427, = NEE 
en ) A I, 
(5.24) 


Equations (5.23) and (5.24) satisfy the conditions 
of (5.17) and (5.18). As a result, a separable filter 
results in separable filters in digital domain and they are 
stable. 


NON Separable tractors 


The first non-separable factor in Chapter IV is 


Pre with a —- -—-l, po dm RT. l as given below 
D 1+30,Ju 
AZ 
or 
1 
H(s} 7S3) Trss, (52255) 
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Fig. 5.3b. The frequency response of 


A? (ez, 71) (ez, 
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and by the double z-transform, we get 


2 -1 -1 
A (1 +Z} ) (1 +z, ) 


zz.) = (5420) 
ee 2 2 2 -1 -1 2 == : 
(^ *4) (4 -4) (z +z, )+ (4 +4)z, Z, 
In the analog domain (5.25b) does not satisfy the 


condition (5.17). We have a singularity at 


wwz = l 

For the second condition, we have a nonessential 
singularity of the second kind at Z| = cd Z, = -]. This 
does not tell us whether it is stable or not, When we 
check the frequency response, it is seen that this factor 
is unstable. 


The second non-separable canonic factor can be 


derived from (4.d) with a q = -1 and Dear l as 





follows. 
ps jos) = rn (52274) 
a o jui 
l + - 
Ju, 
Or 
S 
l 2 
clc E a2 
H(s,1S)) A = 2 F (52,275) 
l + — 
x 
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Then, (5.27b) becomes after performing the double z-transform, 


as follows: 


Q 217) Q7 221) 
ar 
ME 


Tz Kar) 





1122) = 
2(l-z 


This canonic factor does not satisfy either 
stability condition. Hence it is unstable. 
The other series of canonic factors can be derived 


from (4.e). 


For example, with a = -l, ee l, we 
get 
ANS = un (5.294) 
Jj 2 1+ju, +ju, 
or 
Die „os. =P u (5.29b) 
T 2 l +s, +s, : 


or in the digital domain 


(217) (1 zl) 


T(2,,2,) = I— (5.30) 


| NL M x -l- 
(4*4) * A (Z4 +2, ) + (A-4) 2, z 


This canonic factor satisfies the first condition but 
it has a nonessential singularity of the second kind, at 


Z, = 2, = -l and so we should check the frequency 
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response. We see that there is a singularity at 
(07/07) = (1,1). Thus this non-separable canonic factor 
is unstable. 
4. Conclusions 

All of the canonic factors discussed in this 
chapter are collected in Table I. By comparing the results 
shown in this table, we can say that separable stable canonic 
factors result in stable separable filters. As seen from 
Table I, the separable transfer functions l, 2 and 3 which 
start with stable separable ones in the two-dimensional 
continuous frequency domain end up as stable separable trans- 
fer functions in the two-dimensional discrete domain. This 
is also true for the mxn separable canonic factor 4. 

Non-separable canonic factors have very different 
characteristics. The unstable non-separable case ends up 
with a non-separable case in the discrete domain, which is 
also unstable. Although we could not prove this result 
generally, we also could not find any counter example to 
mise. The transfer functions of 5, and 6 start in the 
unstable case, and end up with the unstable case in the 
discrete domain. The transfer functions of 7 starts with 
the stable case in the continuous frequency domain, ends up 
with the unstable case in the discrete domain. Thus we 
have to be very careful with the double z-transform since it 
does not yield stable results every time. It may yield 


unstable results even ıf it starts from the stable case. 
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VI. CONCLUSIONS 


A technique for the use of semi-infinite planes to 
approximate the log magnitude versus log frequency charac- 
teristics of two-dimensional filters is presented. This 
technique is applied to quarter plane filters and works 
well for separable transfer functions. Other non- 
separable canonic (basic) transfer functions are also 
studied in terms of their planar approximation. 

The following separable and nonseparable canonic factors 


have been studied. 


Delle on A DE 

(2) (1 + ae + Ju) al + jp ed. 
EMI c MIS 

Don) EE jo 

EMI go c oO 


The properties and stability of these canonic factors 
are investigated first in the continuous variable domain 
and then after the double bilinear z-transform is applied 
to extend them to two-dimensional recursive digital filters. 
These are all summarized in TableI in Chapter V. An 
important result is that although the double bilinear 


z-transform works well for separable factors, it does not 
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always lead to a stable digital realization starting with a 


stable analog transfer function. 

The following are suggested for further research: 

l. Other types of canonic factors need to be considered. 
For example, the general form: 
H (ju ,307) = [l + (jo, *,)P* (ju41,) S (ju 1.) Go 14) 1? 

2. The work should be extended to non-causal (half- 
plane) filters. 

Jj. Transfer functions which have circular symmetry in 
the discrete domain need to be investigated. 

4. Another approach is to consider half-plane 
recursive algorithms, by working backwards from the discrete 


domain (z to continuous domain (S1,85). 


1122) 
5. Stability of canonic factors needs to be investigated 
in detail. 


6. The use of predistortion should also be considered. 
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